Global solutions of smart building-grid energy management models

ABSTRACT

According to an aspect of the invention, there is provided a method for optimizing a cost of electric power generation in a smart site energy management model, including providing a cost function that models a smart building-grid energy system of a plurality of buildings on a site interconnected with electric power grid energy resources and constraints due to a building model, an electric grid model, and a building-grid interface model, where decision variables for each of the building model, the electric grid model, and the building-grid interface model are box-constrained, and minimizing the cost function subject to the building model constraints, the electric grid model constraints, and building-grid interface model constraints.

CROSS REFERENCE TO RELATED UNITED STATES APPLICATIONS

This application claims priority from “On Solution of a Class of SmartBuilding-Grid Energy Management Models”, U.S. Provisional ApplicationNo. 61/612,570 of Motto, et al., filed Mar. 19, 2012, and “On Solutionof a Class of Smart Building-Grid Energy Management Models”, U.S.Provisional Application No. 61/607,795 of Motto, et al., filed Mar. 7,2012, the contents of both of which are herein incorporated by referencein their entireties.

TECHNICAL FIELD

This disclosure is directed to methods for convex relaxation of a classof smart building-grid energy management models using a semidefiniteprogramming approach.

DISCUSSION OF THE RELATED ART

A smart building-grid energy system is roughly defined as aninterconnection of buildings and electric power grid energy resourceswithin a clearly defined boundary, herein below called a smart site. Asmart site is a strategic or self-interested entity that seeks tomaximize some utility function subject to all applicable technical andbudget constraints. The utility function is a measure of energyefficiency. A smart site may purchase energy from or sell energy to itsexternal environment. On one hand, traditional building energymanagement systems have focused on meeting the energy requirements ofone or more buildings assuming that the buildings are connected to astrong or infinite-capacity utility grid. On the other hand, electricgrid management systems have traditionally modeled relatively largesites as single grid nodes with simple lump load models.

Modeling and solving the coordinated self-interested building andelectric grid management models can help improve energy efficiency andsecurity. From a practical perspective, global solutions to a generalclass of mathematical optimizations are of interest. Such optimizationproblems usually feature nonlinear constraints and many local minima. Inmost current works, local solutions are obtained using heuristicprocedures.

From algorithmic perspective, the field of semidefinite programming(SDP) has received much interest in mathematical optimization in thelast decade. SDP has applications in such diverse fields as traditionalconvex constrained optimization, control theory, and combinatorialoptimization. Because SDPs are solvable via interior-point methods, andusually require about the same amount of computational resources aslinear optimization, most applications can usually be efficiently solvedin practice as well as in theory. There is interest in obtaining tightSDP relaxations of the building-grid energy management optimization andpossibly exploiting the sparsity structure of the system.

Optimization in building energy management system has been much studied.The objective is to meet the thermal and/or electricity demand in thefacility while minimizing a cost function that is a measure of, e.g.,electricity cost, fuel cost, energy losses, etc. Reference V. Chandan,et al., “Modeling and Optimization of a Combined Cooling, Heating andPower Plant System,” in Proc. American Control Conference, 2012, thecontents of which are herein incorporated by reference in theirentirety, presents a building energy management model that is a specialcase of a smart site, specifically a given university campus, andobtained best local solutions using a heuristic procedure. Moreover, theabove methods neither seek a global minimum, nor provide a mechanism tomeasure the gap between the found local minimum and a potential globalminimum.

SUMMARY

Exemplary embodiments of the invention as described herein generallyinclude coordinated building and energy management models, and a convexrelaxation of the ensuing nonlinear, nonconvex model using an interiorpoint method that exploits the sparsity of the model structure for whicha solution yields a tight lower bound. A method according to anembodiment of the invention is potentially large-scale and provides ameans to test whether any solution is global, within engineeringtolerance, and the corresponding optimality gap. A coordinated buildingand grid energy management model according to an embodiment of theinvention has increased energy efficiency potential because acoordinated strategy dominates (weakly or strongly) an uncoordinatedstrategy, a well-known result in strategic game theory, and increasedenergy security potential because a model according to an embodiment ofthe invention does not assume an infinite-capacity grid, but ratherincludes grid operational and security constraints. Hence a feasiblesolution to a model according to an embodiment of the invention may notviolate electric grid operation and security requirements.

According to an aspect of the invention, there is provided a method foroptimizing a cost of electric power generation in a smart site energymanagement model, including providing a cost function ζ that models asmart building-grid energy system of a plurality of buildings on a siteinterconnected with electric power grid energy resources and constraintsdue to a building model, an electric grid model, and a building-gridinterface model, where decision variables for each of the buildingmodel, the electric grid model, and the building-grid interface modelare box-constrained; and minimizing the cost function subject to thebuilding model constraints, the electric grid model constraints, andbuilding-grid interface model constraints, where the building modelincludes a plurality of buildings, one or more gas turbines andgenerators as on-site sources for electricity power, electric chillers,pumps and thermal energy storage tanks, where waste heat in exhaust gasfrom the gas turbines is used to generate steam in one or more heatrecovery steam generator (HRSG) units, steam from the HRSG units drivesone or more absorption chillers to generate cooling, one or more steamturbines to drive additional electricity generators, or provides heatingneeds for the plurality of buildings, the electric grid includes one ormore loads and is modeled as an undirected graph G=(N, E), where N and Edenote sets of electric nodes and branches, respectively, where eachbranch has an origin terminal and a destination terminal, each terminalis connected to an electric node, a voltage at node n₁ is a complexnumber defined by its magnitude and angle, and a power flow at terminalt₁ is a complex number defined by its real and imaginary parts, and thebuilding-grid interface model constrains a building cooling load to besatisfied by chiller operations at all times, and an electricity load tobe satisfied by on-site generation and the electric grid.

According to a further aspect of the invention, the cost function is

${Ϛ = {\sum\limits_{k}\;{\sum\limits_{j}\;\left\lbrack {\rho_{jk}^{T_{j}}{u_{djk}^{T_{j}}\left( x_{jk}^{T_{j}} \right)}} \right\rbrack^{2}}}},$where ρ is a vector of parameters of components of the building modeland electric grid model, u_(d) is a basis vector of d-degreepolynomials, where a dimension of ρ and u(x) is given by

${{u_{d}} = \begin{pmatrix}{{x} + d} \\d\end{pmatrix}},$x is a vector of decision variables of the building model and electricgrid model, j is an index over the components of the building model andelectric grid model, T_(j) denotes one of the components of the buildingmodel and electric grid model, and k is an index of time measurements ofthe parameters and decision variable values, and the constraints areP^(BL)(x)=0, P^(EG)(x)=0, P^(IF)(x)=0, and x≧0, where P^(BL)(x),P^(EG)(x), and P^(IF)(x) are sets of polynomials defining to thebuilding (BL) model constraints, the electric grid (EG) modelconstraints, and the building-grid interface model (IF) constraints,respectively.

According to a further aspect of the invention, the components of thebuilding model include the electric and absorption chillers, a group ofchillers, a group of the plurality of buildings, the generators, the gasturbines, the HRSG units, the pumps, the steam turbines, the thermalenergy storages, and the components of the electric grid model includeloads, nodes, and terminals.

According to a further aspect of the invention, the cost function isreformulated in a semi-definite programming format as

${\min\limits_{X}{C \cdot X}},{{s.t.\mspace{14mu}{{??}^{BL}(X)}} = b^{BL}}$??^(EG)(X) = b^(EG) ??^(IF)(X) = b^(IF) X ≽ 0, rank(X) = 1,by reformulating non-linear terms in the cost function and constraintsas bilinear terms and using a matrix inner product operator on thebilinear terms, where C is a matrix formed of the parameters of thecomponents of the building model and electric grid model, X is a matrixformed of the decision variables of the building model and electric gridmodel, A^(BL)(X), A^(EG)(X), and A^(IF)(X) are reformulations of thebuilding (BL) model constraints, the electric grid (EG) constraints, andthe building-grid interface model (IF) constraints, respectively, in asemi-definite programming format, b^(BL), b^(EG), and b^(IF) areconstants derived from the constraints of the building model andelectric grid model, and

$X = {\begin{bmatrix}1 & x^{T} \\x & {xx}^{T}\end{bmatrix}.}$

According to a further aspect of the invention, the method includesrelaxing the condition rank(X)=1, and solving a following primalsemi-definite program

${\min\limits_{X}{C \cdot X}} - {\mu\;{\ln\left( {\det(X)} \right)}}$subject  to  ??^(BL)(X) = b^(BL), ??^(EG)(X) = b^(EG), ??^(IF)(X) = b^(IF), X ≻ 0,for positive real values of μ as μ approaches zero.

According to a further aspect of the invention, solving the primalsemi-definite program includes imposing optimality conditions on theprimal semi-definite program to deriveC−Z−

^(BL)(y ^(BL))−

^(EG)(y ^(EG))−

^(IF)(y ^(IF))=0,

^(BL)(X)−b ^(BL)=0,

^(EG)(X)−b ^(EG)=0,

^(IF)(X)−b ^(IF)=0,XZ−μI=0,where Z=μX⁻¹, X

0, and Ã^(BL)(y^(BL)), Ã^(EG)(y^(EG)) and Ã^(IF)(y^(IF)) are transposesof A^(BL)(y^(BL)), A^(EG)(y^(EG)) and A^(IF)(y^(IF)) where y^(BL),y^(EQ) and y^(IF) are variables in a dual program of the primalsemi-definite program, initializing a solution (X, y, Z) to an initialpoint (X⁰, y⁰, Z⁰) where X⁰

0,Z⁰

0, choosing a search direction (ΔX, Δy, ΔZ), choosing a primal steplength α_(p) and a dual step length α_(d) that satisfy X+α_(p)ΔX

0 and Z+α_(d)ΔZ

0, and updating X←X+α_(p)ΔX and (y, Z)←α_(d)(Δy, ΔZ).

According to a further aspect of the invention, the method includesrepeating the steps of choosing a search direction, choosing a primalstep length and a dual step length, and updating a current iterate (X,y, Z) until the current iterate satisfies a stopping condition.

According to a further aspect of the invention, choosing a searchdirection comprises solvingΔZ+

^(BL)(Δy ^(BL))+

^(EG)(Δy ^(EG))+

^(IF)(Δy ^(IF))=−r ^(X),

^(BL)(ΔX)=−r ^(BL),

^(EG)(ΔX)=−r ^(EG),

^(IF)(ΔX)=−r ^(IF),

_(P)((ΔX)Z+XΔZ)=−r ^(XZ),wherer ^(X) =C−ΔZ−

^(BL)(Δy ^(BL))−

^(EG)(Δy ^(EG))−

^(IF)(Δy ^(IF)),r ^(BL)=

^(BL)(X)−b ^(BL),r ^(EG)=

^(EG)(X)−b ^(EG),r ^(IF)=

^(IF)(X)−b ^(IF),r ^(XZ)=

_(P)(XZ)−σμI,σε[0,1],H _(P)(XZ)=½(PXZP ⁻¹+(PXZP ⁻¹)^(T)), andH _(P)(ΔXZ+XΔZ)=½(P(ΔXZ+XΔZ)P ⁻¹+(P(ΔXZ+XΔZ)P ⁻¹)^(T)) for anon-singular matrix P.

According to another aspect of the invention, there is provided a methodfor optimizing a cost of electric power generation in a smart siteenergy management model, including providing a semidefinite program

${\min\limits_{X}{C \cdot X}} - {\mu\;{\ln\left( {\det(X)} \right)}}$subject  to  ??^(BL)(X) = b^(BL), ??^(EG)(X) = b^(EG), ??^(IF)(X) = b^(IF), X ≻ 0,where C·X is a cost function that models a smart building-grid energysystem of a plurality of buildings on a site interconnected withelectric power grid energy resources where C is a matrix formed ofparameters of components of a building model and an electric grid model,X is a matrix formed of decision variables of the building model andelectric grid model, A^(BL) (X), A^(EG)(X), and A^(IF)(X) areconstraints due to a building (BL) model, an electric grid (EG) model,and an building-grid interface model (IF), respectively, in asemi-definite programming format, b^(BL), b^(EG), and b^(IF) areconstants derived from the constraints of the building model andelectric grid model, and

${X = \begin{bmatrix}1 & x^{T} \\x & {xx}^{T}\end{bmatrix}},$where a rank condition on X has been relaxed, and solving

${{\min\limits_{X}{C \cdot X}} - {\mu\;{\ln\left( {\det(X)} \right)}}},$subject to the constraints, for positive real values of μ as μapproaches zero.

According to a further aspect of the invention, C·X is derived byreformulating non-linear terms in a function C and constraintsP^(BL)(x)=0, P^(EG) (x)=0, P^(IF) (x)=0, and x≧0, as bilinear terms,using a matrix inner product operator on the bilinear terms toreformulate the bilinear terms in a semi-definite programming format,where

${\varsigma = {\sum\limits_{k}{\sum\limits_{j}\left\lbrack {\rho_{jk}^{T_{j}}{u_{djk}^{T_{j}}\left( x_{jk}^{T_{j}} \right)}} \right\rbrack^{2}}}},$where ρ is a vector of parameters of components of the building modeland electric grid model, u_(d) is a basis vector of d-degreepolynomials, where a dimension of ρ and u(x) is given by

${{u_{d}} = \begin{pmatrix}{{x} + d} \\d\end{pmatrix}},$x is a vector of decision variables of the building model and electricgrid model, j is an index over the components of the building model andelectric grid model, T_(j) denotes one of the components of the buildingmodel and electric grid model, and k is an index of time measurements ofthe parameters and decision variable values, and P^(BL)(x), P^(EG)(x),and P^(IF)(x) are sets of polynomials defining the building (BL) modelconstraints, the electric grid (EG) model constraints, and thebuilding-grid interface model (IF) constraints, respectively.

According to a further aspect of the invention, decision variables foreach of the building model, the electric grid model, and thebuilding-grid interface model are box-constrained.

According to a further aspect of the invention, the building modelincludes a plurality of buildings, one or more gas turbines andgenerators as on-site sources for electricity power, electric chillers,pumps and thermal energy storage tanks, where waste heat in exhaust gasfrom the gas turbines is used to generate steam in one or more heatrecovery steam generator (HRSG) units, steam from the HRSG units drivesone or more absorption chillers to generate cooling, steam turbines todrive additional electricity generators, or provides heating needs forthe plurality of buildings, the electric grid includes one or more loadsand is modeled as an undirected graph G=(N, E), where N and E denotesets of electric nodes and branches, respectively, where each branch hasan origin terminal and a destination terminal, each terminal isconnected to an electric node, a voltage at node n₁ is a complex numberdefined by its magnitude and angle, and a power flow at terminal t₁ is acomplex number defined by its real and imaginary parts, and thebuilding-grid interface model constrains a building cooling load to besatisfied by chiller operations at all times, and an electricity load tobe satisfied by on-site generation and the electric grid.

According to another aspect of the invention, there is provided anon-transitory program storage device readable by a computer, tangiblyembodying a program of instructions executed by the computer to performthe method steps for optimizing a cost of electric power generation in asmart site energy management model.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 illustrates a smart site energy management system architecture,according to an embodiment of the invention.

FIG. 2 illustrates a typical Combined Cooling, Heating and Power (CCHP)system, according to an embodiment of the invention.

FIG. 3 illustrates a unified electric grid branch model, according to anembodiment of the invention.

FIGS. 4(a)-(b) illustrates up-to-eighth-order Taylor seriesapproximations of the cos(x) and sin(x), respectively, over thetheoretical range xε[−π, π], according to an embodiment of theinvention.

FIG. 5 is a block diagram of an exemplary computer system forimplementing a method for convex relaxation of an integrated buildingand energy management model, according to an embodiment of theinvention.

FIG. 6 is a flowchart of a method for solving an SDP relaxation of anintegrated building and energy management model, according to anembodiment of the invention.

DETAILED DESCRIPTION OF EXEMPLARY EMBODIMENTS

Exemplary embodiments of the invention as described herein generallyinclude systems and methods for convex relaxation of an integratedbuilding and energy management model. Accordingly, while the inventionis susceptible to various modifications and alternative forms, specificembodiments thereof are shown by way of example in the drawings and willherein be described in detail. It should be understood, however, thatthere is no intent to limit the invention to the particular formsdisclosed, but on the contrary, the invention is to cover allmodifications, equivalents, and alternatives falling within the spiritand scope of the invention.

Notation:

The notation used throughout the present disclosure is summarized in theNOMENCLATURE Appendix following this specification. Most matricesoccurring herein below will be real symmetric matrices of order n:S_(n)denotes the space of such matrices. The notation x·y denotes the innerproduct of two equi-dimension vectors, defined by x^(T) y. The notationM₁·M₂ denotes the inner product between two such matrices, defined bytr(M^(T) ₁ ^(T)M₂). The associated norm is the Frobenius norm, written∥M∥_(F):=(trM^(T)M)^(1/2) or just ∥M∥, while ∥M∥₂ denotes theL₂-operator norm of a matrix. Norms on vectors will always be Euclideanunless otherwise noted.

The notation X

0 (X

0) means that X is positive semidefinite (positive definite). Thematrices involved in these binary operations are always symmetricmatrices. The notation S₊ ^(n) (S₊₊ ^(n)) denotes the set of positivesemidefinite (positive definite) symmetric matrices of order n. Thenotation X

Y or Y

X means that Y−X

0. The notation X

Y or Y

X means that Y−X

0. If X

0, then write X^(1/2) for the square root of X.

The notation diag(X) denotes the vector of diagonal entries of X. Xdenotes the diagonal matrix with the vector x as its diagonal. Thisnotation can be extended to general block diagonal matrices: if M₁, M₂,. . . , M_(k)εS^(n), then diag(M₁, . . . , M_(k)) denotes the blockdiagonal matrix with the M_(i)'s on its diagonal. Lower-case Romanletters denote vectors and uppercase letters denote matrices.

Additional symbols used here are dual variables, which will bemnemonically displayed on top of the corresponding equality orinequality symbols. If x is a variable, then x and x denote its lowerlimit and upper limit, respectively. Other symbols with a narrower scopeare defined in the subsections where they are used. The symbol π_(ijk)^(CH) denotes the i-th parameter of the j-th chiller (i.e. CH) at timestep k; the symbol x_(ijk) ^(CH) denotes the i-th decision of the j-thchiller (i.e. CH) at time step k; and similarly for the otherbuilding-electric-grid resources; e.g. pump (PM); generator (GN), andnode (ND). In particular, the letter k is a time index, and a parameteror variable subscripted by k means a measurement of the correspondingvalue at time k, and sums over k are sums over measurements of thecorresponding subscripted term at times k.

I. Optimization Model Formulation

FIG. 1 depicts an exemplary, non-limiting smart site energy managementsystem architecture, according to an embodiment of the invention. Asmart site is roughly defined as an interconnection of buildings andelectric power grid energy resources within a clearly defined boundary.A smart site is a strategic or self-interested entity that seeks tomaximize some utility function subject to all applicable technical andbudget constraints. The utility function is a measure of energyefficiency. The smart site 10 according to an embodiment of theinvention depicted in FIG. 1 includes human-machine interfaces (HMIs) 11that visualize results, data analytics engines and smart applicationservers 12, distributed communications 13, and connections 14 to theutility energy grid via auto-reclosers.

I.1. Building Model Constraints

Embodiments of the invention consider a small Combined Cooling, Heatingand Power (CCHP) system that provides heating, cooling and electricitypower to a single building or complex that has a group of severalbuildings. CCHP is a popular technology that integrates cooling, heatingand power generation capabilities on one site. The keystone of CCHP isthat waste heat from a prime mover, such as a gas turbine, can beutilized to meet a variety of thermal and power needs in the complex.

A schematic plot of an exemplary, non-limiting CCHP system is shown inFIG. 2. A gas turbine (GT) and generator (G) serve as the primaryon-site source for electricity power. The waste heat in the exhaust gasfrom gas turbine is used to generate steam in heat recovery steamgenerator (HRSG) unit. Steam from the HRSG unit can either drive theabsorption chiller (AC) to generate cooling, or drive a steam turbine(ST) to drive an additional electricity generator (G), or provideheating needs in the facility. Electric chillers (EC) and thermal energystorage (TES) tanks are also typical in such systems.

Optimization of a CCHP system is often performed on a steady statemodel, assuming that transients die out very fast. Chandan, et al,referenced above, proposed a collection of reduced order models forcomponents such as chillers, TES, HRSG, gas turbine for the purpose ofoptimization.

However, embodiments of the invention can eliminate all task specificfeatures and retain only the functional structures from the models inChandan, et al. to abstract the following set of constraints on modelvariables. These constraints are used to formulate an optimization for ageneral class of building/CCHP systems. In the model equations below,the notation

k refers to a k^(th)-order polynomial whose coefficients can bedetermined from a regression analysis. Note that these model equationsare exemplary and non-limiting, and other embodiments of the inventioncan include other components models and other sets of constraints.

(1) Chiller (CH) Model: These equations describe a typical modelaccording to an embodiment of the invention for chillers in a variableprimary flow system. Under certain simplifying assumptions, the amountof cooling produced by a chiller is determined from chilled water massflow rate through the chiller, expressed in EQ. (1a). The correspondingelectricity consumption is determined by EQ. (1b).x _(3ck) ^(CH) =x _(4ck) ^(CH)π_(3wk) ^(CH)(x _(1gk) ^(CH) −x _(2gk)^(CH))  1ax _(1ck) ^(CH)=

₁(x _(3ck) ^(CH))  (1b)

(2) Thermal Storage (TS) Model: The following constraints describe astratified, two-layered model according to an embodiment of theinvention for thermal energy storage tank. The constraints of EQS. (2)model the charging mode requirements. The constraints of EQS. (3)enforce the discharge mode requirements. The thermodynamics of the twolayers are described by a pair of ordinary differential equations;energy and mass conservation at mixing valves are included as equalityconstraints.

$\begin{matrix}{\mspace{79mu}{x_{7{tk}}^{TS} = {x_{2{gk}}^{CG} - x_{1{mk}}^{CM}}}} & \left( {2a} \right) \\{{\pi_{1{wk}}^{WR}\pi_{3{wk}}^{WR}\frac{\mathbb{d}x_{1{tk}}^{TS}}{\mathbb{d}t}} = {{\pi_{3{tk}}^{TS}x_{7{tk}}^{TS}{\pi_{3{wk}}^{WR}\left( {x_{2{tk}}^{TS} - x_{1{tk}}^{TS}} \right)}} + {\pi_{7{tk}}^{TS}{\pi_{9{tk}}^{TS}\left( {x_{2{tk}}^{TS} - x_{1{tk}}^{TS}} \right)}}}} & \left( {2b} \right) \\{{\pi_{1{wk}}^{WR}\pi_{3{wk}}^{WR}\frac{\mathbb{d}x_{2{tk}}^{TS}}{\mathbb{d}t}} = {{\pi_{8{tk}}^{TS}x_{7{tk}}^{TS}{\pi_{3{wk}}^{WR}\left( {x_{3{tk}}^{TS} - x_{2{tk}}^{TS}} \right)}} + {\pi_{7{tk}}^{TS}{\pi_{9{tk}}^{TS}\left( {x_{1{tk}}^{TS} - x_{2{tk}}^{TS}} \right)}}}} & \left( {2c} \right) \\{\mspace{79mu}{x_{3{tk}}^{TS} = {x_{2{mk}}^{CM} - \pi_{1{gk}}^{CG}}}} & \left( {2d} \right) \\{\mspace{79mu}{{x_{2{gk}}^{CG}x_{1{gk}}^{CG}} = {{x_{7{tk}}^{TS}x_{4{tk}}^{TS}} + {x_{1{mk}}^{CM}x_{3{mk}}^{CM}}}}} & \left( {2e} \right) \\{\mspace{79mu}{x_{7{tk}}^{TS} = {x_{1{mk}}^{CM} - x_{2{gk}}^{CG}}}} & \left( {3a} \right) \\{{\pi_{1{wk}}^{WR}\pi_{3{wk}}^{WR}\frac{\mathbb{d}x_{1{tk}}^{TS}}{\mathbb{d}t}} = {{\pi_{4{tk}}^{TS}x_{7{tk}}^{TS}{\pi_{3{wk}}^{WR}\left( {x_{6{tk}}^{TS} - x_{1{tk}}^{TS}} \right)}} + {\pi_{5{tk}}^{TS}{\pi_{9{tk}}^{TS}\left( {x_{2{tk}}^{TS} - x_{1{tk}}^{TS}} \right)}}}} & \left( {3b} \right) \\{{\pi_{1{wk}}^{WR}\pi_{3{wk}}^{WR}\frac{\mathbb{d}x_{2{tk}}^{TS}}{\mathbb{d}t}} = {{\pi_{6{tk}}^{TS}x_{7{tk}}^{TS}{\pi_{3{wk}}^{WR}\left( {x_{1{tk}}^{TS} - x_{2{tk}}^{TS}} \right)}} + {\pi_{5{tk}}^{TS}{\pi_{9{tk}}^{TS}\left( {x_{1{tk}}^{TS} - x_{2{tk}}^{TS}} \right)}}}} & \left( {3c} \right) \\{\mspace{79mu}{{x_{1{mk}}^{CM}x_{2{mk}}^{CM}} = {{x_{7{tk}}^{TS}x_{5{tk}}^{TS}} + {x_{2{gk}}^{CG}x_{1{gk}}^{CG}}}}} & \left( {3d} \right) \\{\mspace{79mu}{x_{6{tk}}^{TS} = {x_{1{gk}}^{CG} = x_{3{mk}}^{CM}}}} & \left( {3e} \right)\end{matrix}$

(3) Gas Turbine (GT) Model: The constraints of EQS. (4) are data-drivenmodels according to embodiments of the invention for selected input andoutputs of the gas turbine, obtained through regression analysis. Theinput variable is the desired electrical power produced by the gasturbine. The output variables of interest are the corresponding valuesof natural gas mass flow rate, exhaust gas mass flow rate and turbineexit temperature.x _(7tk) ^(TS)=

₁(x _(3gk) ^(GT))  (4a)x _(1gk) ^(GT)=

₂(x _(3gk) ^(GT))  (4b)x _(2gk) ^(GT)=

₄(x _(3gk) ^(GT))  (4c)

(4) Heat Recovery Steam Generator (HR): The HRSG includes multiplecomponents such as super-heater, economizer, boilers, etc. Forsimplicity of exposition, embodiments of the invention can model theHRSG as a lumped, counter-flow heat exchanger.

$\begin{matrix}{x_{3{hk}}^{HR} = {{\mathbb{P}}_{1}\left( x_{6{hk}}^{HR} \right)}} & \left( {5a} \right) \\{\pi_{2{wk}}^{WR} = {{\mathbb{P}}_{1}\left( x_{6{hk}}^{HR} \right)}} & \left( {5b} \right) \\{x_{5{hk}}^{HR} = {\pi_{1{hk}}^{HR}x_{7{hk}}^{HR}}} & \left( {5c} \right) \\{{x_{2{gk}}^{GT} - x_{4{hk}}^{HR}} = {\pi_{2{hk}}^{HR}\left( {x_{2{gk}}^{GT} - x_{1{hk}}^{HR}} \right)}} & \left( {5d} \right) \\{{{\overset{\sim}{x}}_{3{gk}}^{HR}x_{8{hk}}^{HR}} = {x_{1{gk}}^{GT}\left( {x_{2{gk}}^{GT} - x_{4{hk}}^{HR}} \right)}} & \left( {5e} \right) \\{x_{2{hk}}^{HR} = \left\{ \begin{matrix}{0,} & {{\overset{\sim}{x}}_{3{gk}}^{HR} < {\overset{\sim}{x}}_{1{gk}}^{HR}} \\{{x_{3{hk}}^{HR} + {\frac{1}{\pi_{1{gk}}^{GS}}\left( {{\overset{\sim}{x}}_{3{gk}}^{HR} - {\overset{\sim}{x}}_{2{gk}}^{HR}} \right)}},} & {{\overset{\sim}{x}}_{3{gk}}^{HR} \geq {\overset{\sim}{x}}_{2{gk}}^{HR}} \\{x_{3{hk}}^{HR},} & {o.w.}\end{matrix} \right.} & \left( {5f} \right) \\{{\overset{\sim}{x}}_{1{gk}}^{HR} = {\pi_{3{wk}}^{WR}\left( {x_{1{hk}}^{HR} - x_{3{hk}}^{HR}} \right)}} & \left( {5g} \right) \\{{\overset{\sim}{x}}_{2{gk}}^{HR} = {{\pi_{3{wk}}^{WR}\left( {x_{1{hk}}^{HR} - x_{3{hk}}^{HR}} \right)} + \pi_{2{wk}}^{WR}}} & \left( {5h} \right)\end{matrix}$

(5) Steam Loop: The steam loop may be defined as the section of the CCHPsystem associated with the co-generated steam used to meet heatingdemands and drive steam turbines for power production. A simplifiedmodel according to an embodiment of the invention for the steam loop isgiven by EQS. (6), where x_(8hk) ^(HR) denotes the mass flow rate (kg/s)of steam used to drive steam turbine, and x_(1sk) ^(ST) denotes theelectrical power produced (MW) by steam turbine.

$\begin{matrix}{x_{2{sk}}^{ST} = {x_{3{sk}}^{ST}x_{8{hk}}^{HR}}} & \left( {6a} \right) \\{x_{4{mk}}^{CM} = {\left( {1 - x_{3{sk}}^{ST}} \right)x_{8{hk}}^{HR}}} & \left( {6b} \right) \\{{\pi_{1{wk}}^{WR}x_{11k}^{PM}} = {x_{2{sk}}^{ST}\left( {\pi_{1{dk}}^{DA} - \pi_{1{ck}}^{CD}} \right)}} & \left( {6c} \right) \\{x_{12k}^{PM} = {\frac{x_{8{hk}}^{HR}}{\pi_{1{wk}}^{WR}}\left( {x_{7{hk}}^{HR} - \pi_{1{dk}}^{DA}} \right)}} & \left( {6d} \right) \\{x_{1{sk}}^{ST} = {10^{- 3}x_{2{sk}}^{ST}\pi_{1{gk}}^{GS}\pi_{1{sk}}^{ST}{x_{2{hk}}^{HR}\left( {1 - \left( {\pi_{1{ck}}^{CD}/x_{5{hk}}^{HR}} \right)^{\frac{1}{4}}} \right)}}} & \left( {6e} \right)\end{matrix}$

(6) Variable Constraints: All building model decision variables inaccording to embodiments of the invention are box constrained, in thatthey have lower and upper limits.x _(igk) ^(CG) ≦x _(igk) ^(CG) ≦ x _(igk) ^(CG) ,∀i,∀g,∀k  (7a)x _(ick) ^(CH) ≦x _(ick) ^(CH) ≦ x _(ick) ^(CH) ,∀i,∀c,∀k  (7b)x _(imk) ^(CM) ≦x _(imk) ^(CM) ≦ x _(imk) ^(CM) ,∀i,∀m,∀k  (7c)x _(igk) ^(GT) ≦x _(igk) ^(GT) ≦ x _(igk) ^(GT) ,∀i,∀g,∀k  (7d)x _(ihk) ^(HR) ≦x _(ihk) ^(HR) ≦ x _(ihk) ^(HR) ,∀i,∀h,∀k  (7e)x _(ipk) ^(PM) ≦x _(ipk) ^(PM) ≦ x _(ipk) ^(PM) ,∀i,∀p,∀k  (7f)x _(isk) ^(ST) ≦x _(isk) ^(ST) ≦ x _(isk) ^(ST) ,∀i,∀s,∀k  (7g)x _(itk) ^(TS) ≦x _(itk) ^(TS) ≦ x _(itk) ^(TS) ,∀i,∀t,∀k  (7h)I.2. Electric Grid Model Constraints

An electric grid (EG) can be modeled as a standard undirected graphG=(N, E), where N and E denote the sets of electric nodes and branches,respectively. A branch has two terminals: origin and destinationterminals. Each terminal is connected to an electric node. Examples ofelectric branches are underground or overhead lines, on-load taptransformers, phase-shifting transformers, series reactors and seriescapacitors. In concept, generators, loads, and shunt devices, such asshunt capacitors and reactors, can be modeled as branch devices with oneof their associated terminals connected to ground.

FIG. 3 depicts a unified branch model of an electric grid according toan embodiment of the invention, whereby the subsequent equations can bereadily derived. On the left, a generator s_(k) is connected to node k,and load s₁ is connected to node l on the right. The branch modelincludes transformer taps between nodes k and k′, and between nodes l′and l, with branches connecting to internal nodes k′ and l′. The voltageat node n₁ is a complex number defined by its magnitude and angle, heredenoted x_(1n) ₁ ^(ND) and x_(2n) ₁ ^(ND), respectively. The power flowat terminal t₁ is a complex number defined by its real and imaginaryparts, here denoted x_(1n) ₁ ^(TE) and x_(2n) ₁ ^(TE), respectively. Thesymbol π₁ ^(BR) denotes the series admittance of the branch connected toterminals t₁ and t₂; π₂ ^(BR) and π₃ ^(BR) model the shunt admittancesof the branch, and π₄ ^(BR) and π₅ ^(BR) model the transformation ratiobetween the branch external and internal nodes, at the origin anddestination nodes, respectively.

The aforementioned grid resource model according to an embodiment of theinvention is sufficiently general to include practical management modelsof some distributed energy resources such as PV generators and theirinverter interfaces, as well as batteries. For example, several PV andbattery resources management models known in the art are linear and aresupported by a general framework according to an embodiment of theinvention.

Terminal Constraints:

According to an embodiment of the invention, the following set ofconstraints are defined for all terminals t, at all times k. The flow ata terminal t is defined in terms of the voltage of the node associatedwith the terminal, here denoted n(t), as well as the voltage of the nodeassociated with a terminal, referred to herein as its branch-imageterminal, here denoted n({tilde over (t)}).x _(1tk) ^(TE)=π_(1tk) ^(TE) x _(1n(t)k) ^(ND) x _(1n(t)k) ^(ND) +x_(1n(t)k) ^(ND) x _(1n({tilde over (t)})k) ^(ND){π_(2tk) ^(TE) cos(x_(1n(t)k) ^(ND) −x _(1n({tilde over (t)})k) ^(ND)+π_(4tk) ^(TE) sin(x_(1n(t)k) ^(ND) −x _(1n({tilde over (t)})k) ^(ND))}  (8a)x _(1tk) ^(TE)=π_(1tk) ^(TE) x _(1n(t)k) ^(ND) x _(1n(t)k) ^(ND) +x_(1n(t)k) ^(ND) x _(1n({tilde over (t)})k) ^(ND){π_(2tk) ^(TE) cos(x_(1n(t)k) ^(ND) −x _(1n({tilde over (t)})k) ^(ND))+π_(4tk) ^(TE) sin(x_(1n(t)k) ^(ND) −x _(1n({tilde over (t)})k) ^(ND))}  (8a)x _(2tk) ^(TE)=−π_(3tk) ^(TE) x _(1n(t)k) ^(ND) x _(1n(t)k) ^(ND) +x_(1n(t)k) ^(ND) x _(1n({tilde over (t)})k) ^(ND){π_(2tk) ^(TE) sin(x_(1n(t)k) ^(ND) −x _(1n({tilde over (t)})k) ^(ND))−π_(4tk) ^(TE) cos(x_(1n(t)k) ^(ND) −x _(1n({tilde over (t)})k) ^(ND))}  (8b)(x _(1tk) ^(TE))²+(x _(2tk) ^(TE))²≦π_(1tk) ^(TE)  (8c)

The constraints of EQ. (8a) define the real part of the power flow atterminal t at time k. The constraints of EQ. (8b) define the imaginarypart of the power flow at terminal t at time k. The constraints of EQ.(8c) enforce the thermal rating at terminal t at time k. Theseconstraints contain nonlinear terms involving the sine and cosinefunctions. Here, embodiments of the invention replace these functions bypolynomial approximations using the Taylor series expansions, up toseventh order. FIGS. 4(a)-(b) in illustrates up-to-eighth-order Taylorseries approximations of the cos(x) and sin(x), respectively, over thetheoretical range xε[−π,π], which largely covers the angle range for allpractical operational and planning purposes, according to an embodimentof the invention. FIGS. 4(a)-(b) shows the quality of thisapproximation, respectively, for the sine (P₁, P₃, P₅, and P₇) andcosine (P₂, P₄, P₆, and P₈), the seventh order of each of which ishighly-accurate for most realistic grid operation and planning modelingpurposes. According to embodiments, the node angle differences are morerestricted so that lower order polynomial approximations may beappropriate, thereby decreasing the size of the optimization model.

The constraints of EQS. (8a) and (8b) introduce hard nonlinearity andnon-convexity into the feasible region of grid management modelsaccording to embodiments of the invention. However, note that theseconstraints exhibit sparsity, owing to the dependency of x_(tk) ^(TE) ontwo 2-dimensional variables, x_(n(t)) ^(ND) and x_(n({tilde over (t)}))^(ND).

Node Constraints:

According to embodiments, the following set of constraints are definedfor all nodes n, at all times k.

$\begin{matrix}{{{\sum\limits_{\{{{g|{n{(g)}}} = n}\}}x_{1{gk}}^{GN}} - {\sum\limits_{\{{{l|{n{(l)}}} = n}\}}x_{1{lk}}^{LD}} - {\sum\limits_{{t|{n{(t)}}} = n}x_{1{tk}}^{TE}}} = 0} & \left( {9a} \right) \\{{{\sum\limits_{\{{{g|{n{(g)}}} = n}\}}x_{2{gk}}^{GN}} - {\sum\limits_{\{{{l|{n{(l)}}} = l}\}}x_{2{lk}}^{LD}} - {\sum\limits_{{t|{n{(t)}}} = n}x_{2{tk}}^{TE}}} = 0} & \left( {9b} \right)\end{matrix}$where n(t) refers to all nodes n associated with a particular terminalt. The constraints of EQ. (9a) state that, for every node n, at everytime k, the real part of the total power produced minus the real part ofthe total power consumed is equal to the net power injected to theelectric grid. The constraints of EQ. (9b) state that, for every node n,at every time k, the imaginary part of the total power produced minusthe imaginary part of the total power consumed is equal to the net powerinjected to the electric grid.

Note that EQS. (9) are highly sparse since the degree of nodes inreal-world electric grids is typically small (1 to 5 branches per node).

Variable Constraints:

According to embodiments of the invention, all electric grid decisionvariables are box constrained; that is, they admit lower and upperbounds.x _(igk) ^(GN) ≦x _(igk) ^(GN) ≦ x _(igk) ^(GN) ,∀i,∀g,∀k  (10a)x _(ilk) ^(LD) ≦x _(ilk) ^(LD) ≦ x _(ilk) ^(LD) ,∀i,∀l,∀k  (10b)x _(itk) ^(TE) ≦x _(itk) ^(TE) ≦ x _(itk) ^(TE) ,∀i,∀t,∀k  (10c)x _(ink) ^(ND) ≦x _(ink) ^(ND) ≦ x _(ink) ^(ND) ,∀i,∀n,∀k  (10d)I.3. Building-Grid Interface Constraints

Embodiments of the invention express the constraints that link buildingand electric grid management models as follows:

$\begin{matrix}{\mspace{79mu}{{{\sum\limits_{c}x_{3{ck}}^{CH}} = {\sum\limits_{m}\pi_{3{mk}}^{CM}}},{\forall k}}} & \left( {11a} \right) \\{{{{\sum\limits_{m}\pi_{5{mk}}^{CM}} + {\sum\limits_{g}\pi_{1{gk}}^{GT}} + {\sum\limits_{s}x_{1{sk}}^{ST}}} = {{\sum\limits_{m}\pi_{2{mk}}^{CM}} + {\sum\limits_{p}x_{1{pk}}^{PM}} + {\sum\limits_{c}x_{1{ck}}^{CH}}}},\mspace{20mu}{\forall k}} & \left( {11b} \right)\end{matrix}$EQ. (11a) requires that the campus cooling load is satisfied by theoperation of chillers at all time. EQ. (11b) ensures that theelectricity load on campus is satisfied by on-site generation and powergrid.I.4. Objective Function

Embodiments of the invention may assume a generic objective function tobe minimized, which may be cast as a piecewise linear function or a sumof squares (SOS) polynomial. A polynomial pεR_(2d)[x] is an SOSpolynomial if there exists k polynomials, {{tilde over (p)}_(i)}_(i=1)^(k)εR[x], such that p(x)=Σ_(i=1) ^(k){tilde over (p)}_(i)(x)². LetS_(2d)[x] denote the set of SOS polynomial in x of degree at most 2d.For example, p(x₁,x₂)=(x₁ ²−3x₂−1)²+(2x₁x₂+x₁−5)² is an SOS polynomialof degree 4, with k=2, {tilde over (p)}₁(x₁,x₂)=x₁ ²−3x2−1, and {tildeover (p)}₂(x₁,x₂)=2x₁x₂+x₁−5.

Note that a d-degree polynomial pεR_(d)[x] may be cast asp(x)=ρ^(T)u_(d)(x), where ρεR^(∥u) ^(d) ^(|) is a vector of parametersπ_(ijk) ^(XY), and u_(d)(x)εR^(∥u) ^(d) ^(|) is a basis vector ofr-degree polynomials. The dimension of ρ and u(x) is given by

${{u_{d}} = \begin{pmatrix}{{x} + d} \\d\end{pmatrix}},$i.e., the number of combinations of d in |x|+d. The vector u_(d)(x) maybe expanded, in terms of the basis elements, as u_(d)(x)=vec{1, x₁, x₂,. . . , x_(n), x₁ ², x₁, x₂, x₃, . . . , x₁ ^(r), . . . , x_(n) ^(r)}.For example, if xεR³ and pεP₂[x], one has u₂(x₁, x₂, x₃)=vec{1, x₁, x₂,x₃, x₁ ², x₁x₂, x₁x₃, x₂ ², x₂x₃, x₃ ²}.

Using the aforementioned notation, a set of 2d-degree SOS polynomialscan be formally written as follows:

$\begin{matrix}{{{??}_{2d}\lbrack x\rbrack} = {\left\{ {{{\sum\limits_{j = 1}^{k}\left( {\rho_{j}^{T}{u_{d}(x)}} \right)}:{k \geq 1}},{\rho_{j} \in {\mathbb{R}}^{(\begin{matrix}{{x} + d} \\d\end{matrix})}}} \right\}.}} & (12)\end{matrix}$

A smart site optimization model objective function according to anembodiment of the invention includes decomposable terms of subsets ofbuilding and/or electric grid model variables. Using the aforementionednotation, a smart site objective function, ξ, according to an embodimentof the invention, takes the form

$\varsigma = {\sum\limits_{k}{\sum\limits_{j}\left\lbrack {\rho_{jk}^{T_{j}}{u_{djk}^{T_{j}}\left( x_{jk}^{T_{j}} \right)}} \right\rbrack^{2}}}$where T one of the model components and j represents the correspondingsubscripts, respectively, and may be cast as follows:

$\begin{matrix}{{\zeta = {\sum\limits_{k}\left\{ {{\sum\limits_{g}\left\lbrack {\rho_{gk}^{CG} \cdot {u_{dgk}^{CG}\left( x_{gk}^{CG} \right)}} \right\rbrack^{2}} + {\sum\limits_{c}\left\lbrack {\rho_{ck}^{CH} \cdot {u_{dck}^{CH}\left( x_{ck}^{CH} \right)}} \right\rbrack^{2}} + {\sum\limits_{m}\left\lbrack {\rho_{mk}^{CM} \cdot {u_{dmk}^{CM}\left( x_{mk}^{CM} \right)}} \right\rbrack^{2}} + {\sum\limits_{g}\left\lbrack {\rho_{gk}^{GN} \cdot {u_{dgk}^{GN}\left( x_{gk}^{GN} \right)}} \right\rbrack^{2}} + {\sum\limits_{gl}\left\lbrack {\rho_{gk}^{GT} \cdot {u_{dgk}^{GT}\left( x_{gk}^{GT} \right)}} \right\rbrack^{2}} + {\sum\limits_{h}\left\lbrack {\rho_{hk}^{HR} \cdot {u_{dhk}^{HR}\left( x_{hk}^{HR} \right)}} \right\rbrack^{2}} + {\sum\limits_{l}\left\lbrack {\rho_{lk}^{LD} \cdot {u_{dlk}^{LD}\left( x_{lk}^{LD} \right)}} \right\rbrack^{2}} + {\sum\limits_{n}\left\lbrack {\rho_{nk}^{ND} \cdot {u_{dnk}^{ND}\left( x_{nk}^{ND} \right)}} \right\rbrack^{2}} + {\sum\limits_{p}\left\lbrack {\rho_{pk}^{PM} \cdot {u_{dpk}^{PM}\left( x_{p\; k}^{PM} \right)}} \right\rbrack^{2}} + {\sum\limits_{s}\left\lbrack {\rho_{sk}^{ST} \cdot {u_{dsk}^{ST}\left( x_{sk}^{ST} \right)}} \right\rbrack^{2}} + {\sum\limits_{t}\left\lbrack {\rho_{tk}^{TE} \cdot {u_{dtk}^{TE}\left( x_{tk}^{TE} \right)}} \right\rbrack^{2}} + {\sum\limits_{t}\left\lbrack {\rho_{tk}^{TS} \cdot {u_{dtk}^{TS}\left( x_{tk}^{TS} \right)}} \right\rbrack^{2}}} \right\}}},} & (13)\end{matrix}$where the symbols U_(dgk) ^(CG)(x_(gk) ^(CG)), . . . , u_(dtk)^(TS)(x_(tk) ^(TS)) denote basis vectors of at most d-degree polynomialswith appropriate dimensions, as defined above, e.g., in EQ. (12).

According to embodiments, only a small subset of variables typicallyappear in the objective function. Hence, the objective functionparameters (i.e. ρ_(j)'s) are highly sparse. For example, the cost ofelectricity purchased from the external electric grid can be minimizedtogether with the cost of fuel used to power the gas turbine. For thisexample, the objective function defined in EQ. (13) simplifies asfollows:

$\begin{matrix}{\varsigma = {\sum\limits_{k}\left\{ {{\sum\limits_{m}{\rho_{mk}^{CM}x_{mk}^{CM}}} + {\sum\limits_{g}{\rho_{gk}^{GT}x_{gk}^{GT}}}} \right\}}} & (14)\end{matrix}$

In other energy management applications, the cost of generation isexpressed as a quadratic function in the power produced by thegenerators. It is well-known that any polynomial of degree 2 inxεR^(|x|), |x|≧1 can be cast as a sum of squares polynomials. In thiscase, the objective function is a quadratic:

$\begin{matrix}{\varsigma = {\sum\limits_{k}{\sum\limits_{m}{\left\lbrack {\rho_{gk}^{GN} \cdot {u_{dgk}^{GN}\left( x_{gk}^{GN} \right)}} \right\rbrack^{2}.}}}} & (15)\end{matrix}$I.5 Compact Model

A smart side energy management optimization model according to anembodiment of the invention can minimize the objective function of EQ.(13) subject to the constraints of EQS. (1)-(11). Embodiments of theinvention can recast the aforementioned model using matrix vectornotation. It is conceivable that the matrices are highly sparse formedium- to large-scale real-life smart site models. First, the variousdecision variables may be packed into real vectors of appropriatedimensions:x ^(CG) =vec{x _(1gk) ^(CG) ,x _(2gk) ^(CG)}_(g,k),  (16a)x ^(CH) =vec{x _(1ck) ^(CH) ,x _(2ck) ^(CH) ,x _(3ck) ^(CH) ,x _(4ck)^(CH)}_(c,k),  (16b)x ^(CM) =vec{x _(1mk) ^(CM) ,x _(2mk) ^(CM) ,x _(3mk) ^(CM) ,x _(4mk)^(CM) ,x _(5mk) ^(CM)}_(m,k),  (16c)x ^(GN) =vec{x _(1gk) ^(GN) ,x _(2gk) ^(CN)}_(g,k),  (16d)x ^(GT) =vec{x _(1gk) ^(GT) ,x _(2gk) ^(GT) ,x _(3gk) ^(GT) ,x _(4gk)^(GT)}_(g,k),  (16e)x ^(HR) =vec{x _(1hk) ^(HR) ,x _(2hk) ^(HR) ,x _(3hk) ^(HR) ,x _(4hk)^(HR) ,x _(5hk) ^(HR) ,x _(7hk) ^(HR) ,x _(8hk) ^(HR)}_(h,k),  (16f)x ^(LD) =vec{x _(1lk) ^(LD) ,x _(2lk) ^(LD)}_(l,k),  (16g)x ^(ND) =vec{x _(1nk) ^(ND) ,x _(2nk) ^(ND)}_(n,k),  (16h)x ^(PM) =vec{x _(1pk) ^(PM)}_(p,k),  (16i)x ^(ST) =vec{x _(1sk) ^(ST) ,x _(2sk) ^(ST) ,x _(3sk)^(ST)}_(s,k),  (16j)x ^(TE) =vec{x _(1tk) ^(TE) ,x _(2tk) ^(TE)}_(t,k),  (16k)x ^(TS) =vec{x _(1tk) ^(TS) ,x _(2tk) ^(TS) ,x _(3tk) ^(TS) ,x _(4tk)^(TS) ,x _(5tk) ^(TS) ,x _(6tk) ^(TS) ,x _(7tk) ^(TS) ,x _(8tk) ^(TS) ,x_(9tk) ^(TS)}_(t,k),  (16l)

Finally, embodiments of the invention pack all the smart siteoptimization model decision variables in to a vector x with appropriatedimension:x=vec{x ^(CG) ,x ^(CH) ,x ^(CM) ,x ^(GN) ,x ^(GT) ,x ^(HR) ,x ^(LD) ,x_(n) ^(ND) ,x ^(PM) _(,x) ^(TE) ,x ^(TS)}  (17)

Projection matrices of appropriate dimensions, here denoted P*, may beused to obtain the aforementioned sub-vectors from x. For example,x^(GN)=P^(GN)x, and x^(TE)=P^(TE)x, where P^(GN) and P^(TE) areprojection matrices.

A smart site energy management model according to an embodiment of theinvention can be cast compactly as follows:

$\begin{matrix}{\min\limits_{x}{\varsigma(x)}} & \left( {18a} \right) \\{{s.t.\mspace{14mu}{{??}^{BL}(x)}} = 0} & \left( {18b} \right) \\{{{??}^{EG}(x)} = 0} & \left( {18c} \right) \\{{{??}^{IF}(x)} = 0} & \left( {18d} \right) \\{{\underset{\_}{x} \leq x \leq \overset{\_}{x}},} & \left( {18e} \right)\end{matrix}$where P^(BL)(x), P^(EG)(x), and P^(IF)(x) denote the sets of polynomialsdefining the building (BL) model constraints, the electric grid (EG)model constraints, and the building-grid interface model (IF)constraints, respectively:

^(BL)(x)=vec{p ₁ ^(BL)(x), . . . ,p _(c) _(b) ^(BL)(x)},  (19a)

^(EG)(x)=vec{p ₁ ^(EG)(x), . . . ,p _(c) _(e) ^(EG)(x)},  (19b)

^(IF)(x)=vec{p ₁ ^(IF)(x), . . . ,p _(c) _(i) ^(IF)(x)},  (19c)where c_(b), c_(e), and c_(i) denote the number of building, electricgrid, and interface constraints, respectively. The constraints of EQS(17b), (17c) and (17d) are equivalent to EQS. (1)-(6), (8)-(9), and(11), respectively. The constraints of EQ. (17e) are equivalent to EQS.(7) and (10). The objective function ξ is a linear sum of squarepolynomials in BL and EG decision variables.

Finally, embodiments standardize the model to the following format:

$\begin{matrix}{\min\limits_{x}{\varsigma(x)}} & \left( {20a} \right) \\{{s.t.\mspace{14mu}{{??}^{BL}(x)}} = 0} & \left( {20b} \right) \\{{{??}^{EG}(x)} = 0} & \left( {20c} \right) \\{{{??}^{IF}(x)} = 0} & \left( {20d} \right) \\{x \geq 0.} & \left( {20e} \right)\end{matrix}$which is equivalently derived from EQS. (18) using the followingprocedure:

-   -   (i) change the variable reference: x←x−x, x← x−x;    -   (ii) rewrite the ensuing inequality constraints x≦ x−x as        equality constraints: x+s= x, where s≧0 is an auxiliary        variable;    -   (iii) append s to x; and    -   (iv) append the new equality constraints to EQS. (20b) and        (20c), as appropriate.        II. Semidefinite Relaxation        II.1 Sparse SDP Formulation

A smart site model according to an embodiment of the invention presentedabove is a nonlinear and nonconvex optimization model, possibly with adisconnected feasibility set. The constraints defining the feasibleregion contain linear and nonlinear terms. Embodiments may cast thenonlinear constraints in the smart side model as bilinear terms. Thiscan be achieved by the addition of auxiliary variables wherevernecessary. For example, the constraints of EQ. (4c) are equivalent to:a ₄ x ₂ ² +a ₃ x ₂ x ₁ +a ₂ x ₁ ² +a ₁ x ₁ =b ₁,  (21a)x ₂ =x ₁ ²,  (21b)EQ. (21a) is equivalent to

$\begin{matrix}{{{{\begin{pmatrix}1 \\x_{1} \\x_{2}\end{pmatrix}^{T}\begin{bmatrix}0 & {\frac{1}{2}a_{1}} & 0 \\{\frac{1}{2}a_{1}} & a_{2} & {\frac{1}{2}a_{3}} \\0 & {\frac{1}{2}a_{3}} & a_{4}\end{bmatrix}}\begin{pmatrix}1 \\x_{1} \\x_{2}\end{pmatrix}} = b_{1}},} & \left( {21c} \right)\end{matrix}$or, using the matrix inner product operator, to

$\begin{matrix}{{{\begin{bmatrix}0 & {\frac{1}{2}a_{1}} & 0 \\{\frac{1}{2}a_{1}} & a_{2} & {\frac{1}{2}a_{3}} \\0 & {\frac{1}{2}a_{3}} & a_{4}\end{bmatrix} \cdot \begin{pmatrix}1 \\x_{1} \\x_{2}\end{pmatrix}}\begin{pmatrix}1 \\x_{1} \\x_{2}\end{pmatrix}^{T}} = {b_{1}.}} & \left( {21d} \right)\end{matrix}$

The recast of EQ. (21b) to an equivalent form, using the matrix innerproduct operator, follows from the previous example.

According to embodiments of the invention, applying the scheme in EQ.(21), the smart-site energy management model of EQS. (18) can bereformulated at step 61 in a semi-definite programming format as in EQS.(22):

$\begin{matrix}{{\min\limits_{X}{C \cdot X}},} & \left( {22a} \right) \\{{{s.t.\mspace{14mu}{{??}^{BL}(X)}} = b^{BL}},} & \left( {22b} \right) \\{{{??}^{EG}(X)} = b^{EG}} & \left( {22c} \right) \\{{{??}^{IF}(X)} = b^{IF}} & \left( {22d} \right) \\{{X = \begin{bmatrix}1 & x^{T} \\x & {xx}^{T}\end{bmatrix}},} & \left( {22e} \right)\end{matrix}$or equivalently as EQS. (23):

$\begin{matrix}{\min\limits_{X}{C \cdot X}} & \left( {23a} \right) \\{{s.t.\mspace{14mu}{{??}^{BL}(X)}}\overset{y^{BL}}{=}b^{BL}} & \left( {23b} \right) \\{{{??}^{EG}(X)}\overset{y^{EG}}{=}b^{EG}} & \left( {23c} \right) \\{{{??}^{IF}(X)}\overset{y^{IF}}{=}b^{IF}} & \left( {23d} \right) \\{{X \succcurlyeq 0},} & \left( {23e} \right) \\{{{{rank}(X)} = 1},} & \left( {23f} \right)\end{matrix}$where x denotes a real symmetric matrix aggregating all aforementioneddecision variables, A^(BL)(X), A^(EG)(X), and A^(IF)(X) reformulate thebuilding (BL) model constraints, the electric grid (EG) constraints, andthe building-grid interface model (IF) constraints, respectively, in asemi-definite programming format, and y^(BL), y^(EG), and y^(IF) are therespective dual variables. For example,

^(BL)(X)=vec{A ₁ ^(BL) ·X, . . . ,A _(c) _(b) ^(BL) ·X},  (24a)

^(EG)(X)=vec{A ₁ ^(EG) ·X, . . . ,A _(c) _(e) ^(EG) ·X},  (24b)

^(IF)(X)=vec{A ₁ ^(IF) ·X, . . . ,A _(c) _(i) ^(IF) ·X},  (24c)where c_(b), c_(e), and c_(i) denote the number of building, electricgrid, and interface constraints, respectively. EQS. (22) and (23)represents a change from optimizing over vector spaces to optimizingover the space of symmetric matrices.

Note that EQS. (23) are nonlinear and nonconvex due to the rankcondition on the matrix variable X. Embodiments of the invention mayrelax the rank condition on X which leads to a convex relaxation of asmart site management model according to embodiments of the invention.Embodiments of the invention may solve the ensuing model using aninterior point method that exploits the sparsity structure.

II.2. Solution by Primal-Dual Interior-Point Method

The algebra of a primal-dual interior point method for an SDP relaxationof the management model according to an embodiment of the invention isoutlined as follows. It has been shown that, using a primal-dualinterior point method, SDP can be solved to optimality with an errorbounded by a given tolerance, say ε, in polynomial time.

A relaxed primal SDP according to an embodiment of the invention may bestated as follows:

$\begin{matrix}{\min\limits_{X}{C \cdot X}} & \left( {25a} \right) \\{{{s.t.\mspace{14mu}{{??}^{BL}(X)}}\overset{y^{BL}}{=}b^{BL}},} & \left( {25b} \right) \\{{{{??}^{EG}(X)}\overset{y^{EG}}{=}b^{EG}},} & \left( {25c} \right) \\{{{{??}^{IF}(X)}\overset{y^{IF}}{=}b^{IF}},} & \left( {25d} \right) \\{X \succcurlyeq 0.} & \left( {25e} \right)\end{matrix}$

According to an embodiment of the invention, the dual of EQS. (25) is asfollows:

$\begin{matrix}{{\min\limits_{y,Z}{b^{BL} \cdot y^{BL}}} + {b^{EG} \cdot y^{EG}} + {b^{IF} \cdot y^{IF}}} & \left( {26a} \right) \\{{{{s.t.\mspace{14mu}{{\overset{\sim}{??}}^{BL}\left( y^{BL} \right)}} + {{\overset{\sim}{??}}^{EG}\left( y^{EG} \right)} + {{\overset{\sim}{??}}^{IF}\left( y^{IF} \right)} + Z} = C},} & \left( {26b} \right) \\{Z \succcurlyeq 0.} & \left( {26c} \right)\end{matrix}$where Ã^(BL)(y^(BL)), Ã^(EG)(y^(EG)) and Ã^(IF)(y^(IF)) denotes thetranspose of A^(BL)(y^(BL)), A^(EG)(y^(EG)), and A^(IF)(y^(IF)),respectively:

$\begin{matrix}{{{{\overset{\sim}{??}}^{BL}(y)} = {\sum\limits_{i = 1}^{c_{b}}{y_{i}A_{i}^{BL}}}},} & \left( {27a} \right) \\{{{{\overset{\sim}{??}}^{EG}(y)} = {\sum\limits_{i = 1}^{c_{e}}{y_{i}A_{i}^{EG}}}},} & \left( {27b} \right) \\{{{{\overset{\sim}{??}}^{IF}(y)} = {\sum\limits_{i = 1}^{c_{b}}{y_{i}A_{i}^{IF}}}},} & \left( {27c} \right)\end{matrix}$where c_(b), c_(e), and c_(i) denote the number of building, electricgrid, and interface constraints, respectively.

Given a primal-dual feasible (X, y, Z), i.e. if X and (y, Z) arefeasible solutions for EQS. (25) and (26), respectively, the duality gapis:C·X−b ^(BL) ·y ^(BL) −b ^(EG) ·y ^(EG) −b ^(IF) ·y ^(IF) =X·Z≧0  (28)where equality holds; i.e. Z·X=0, if and only if XZ=0.

A solution approach according to an embodiment of the invention ispresented in FIG. 6 and includes solving the following parameterized SDPproblem for a variety of positive real values μ as μ tends to zero:

$\begin{matrix}{{\min\limits_{X}{C \cdot X}} - {\mu\;{\ln\left( {\det(X)} \right)}}} & \left( {29a} \right) \\{{{s.t.\mspace{14mu}{{??}^{BL}(X)}} = b^{BL}},} & \left( {29b} \right) \\{{{{??}^{EG}(X)} = b^{EG}},} & \left( {29c} \right) \\{{{{??}^{IF}(X)} = b^{IF}},} & \left( {29d} \right) \\{X \succ 0.} & \left( {29e} \right)\end{matrix}$where ln denote the natural logarithm. Referring now to the figure, atstep 60, the optimality or Karush-Kuhn-Tucker (KKT) conditions areimposed on EQS. (29) to yield:C−μX ⁻¹−

^(BL)(y ^(BL))−

^(EG)(y ^(EG))−

^(IF)(y ^(IF))=0,  (30a)

^(BL)(X)−b ^(BL)=0,  (30b)

^(EG)(X)−b ^(EG)=0,  (30c)

^(IF)(X)−b ^(IF)=0,  (30d)X

0.  (30e)

Embodiments of the invention define Z=μX⁻¹, which implies XZ=μI, torecast the KKT conditions as follows:C−Z−

^(BL)(y ^(BL))−

^(EG)(y ^(EG))−

^(IF)(y ^(IF))=0,  (31a)

^(BL)(X)−b ^(BL)=0,  (31b)

^(EG)(X)−b ^(EG)=0,  (31c)

^(IF)(X)−b ^(IF)=0,  (31d)XZ−μI=0,  (31e)where X

0.

A nonlinear algebraic system can be solved using any appropriateNewton-like method, i.e. any method adapted from the following genericprocedure, with reference to the steps shown in FIG. 6:

-   -   1. Select a stopping criteria, and choose an initial point (X⁰,        y⁰, Z⁰), such that X⁰        0 and Z⁰        0. Set (X, y, Z)←(X⁰, y⁰, Z⁰) (Step 61).    -   2. REPEAT    -   3. Choose a search direction (ΔX, Δy, ΔZ) (Step 62).    -   4. Choose a primal step length α_(p) and a dual step length        α_(d) such the following conditions hold: (Step 63)        X+α _(p) ΔX        0.  (32a)        Z+α _(d) ΔZ        0.  (32b)    -   5. Update the variables: (Step 64)        X←X+α _(p) ΔX,  (33a)        (y,Z)←(y,Z)+α_(d)(Δy,Δz).  (33b)    -   6. UNTIL the current iterate (X, y, Z) satisfies the stopping        criteria. (Step 65)

The above steps 1-6 would be repeated from step 66 for decreasing valuesof Z=μX⁻¹.

The search direction (ΔX, Δy, ΔZ) can be determined by solving afollowing algebraic system according to an embodiment of the invention:ΔZ+

^(BL)(Δy ^(BL))+

^(EG)(Δy ^(EG))+

^(IF)(Δy ^(IF))=−r ^(X),  (34a)

^(BL)(ΔX)=−r ^(BL),  (34b)

^(EG)(ΔX)=r ^(EG),  (34c)

^(IF)(ΔX)=−r ^(IF),  (34d)

_(P)((ΔX)Z+XΔZ)=−r ^(XZ),  (34e)wherer ^(X) =C−ΔZ−

^(BL)(Δy ^(BL))−

^(EG)(Δy ^(EG))−

^(IF)(Δy ^(IF)),  (35a)r ^(BL)=

^(BL)(X)−b ^(BL),  (35b)r ^(EG)=

^(EG)(X)−b ^(EG),  (35c)r ^(IF)=

^(IF)(X)−b ^(IF),  (35d)r ^(XZ)=

_(P)(XZ)−σμI  (35e)σε[0,1]  (35f),and H_(P) is a symmetrization operator which, for a non-singular matrixP, ensures that XZ is symmetric at every iterate:H_(P)A=½(PAP⁻¹+(PAP⁻¹)^(T)), thusH _(P)(XZ)=½(PXZP ⁻¹+(PXZP ⁻¹)^(T))  (36a)H _(P)(ΔXZ+XΔZ)=½(P(ΔXZ+XΔZ)P ⁻¹+(P(ΔXZ+XΔZ)P ⁻¹)^(T))  (36b)Exemplary, non-limiting choices for P are P=Z^(1/2) and P=I.

Iterating the above steps until convergence will yield global lowerbound for the cost function.

III. System Implementations

It is to be understood that embodiments of the present invention can beimplemented in various forms of hardware, software, firmware, specialpurpose processes, or a combination thereof. In one embodiment, thepresent invention can be implemented in software as an applicationprogram tangible embodied on a computer readable program storage device.The application program can be uploaded to, and executed by, a machinecomprising any suitable architecture.

FIG. 5 is a block diagram of an exemplary computer system forimplementing a method and system for convex relaxation of an integratedbuilding and energy management model, according to an embodiment of theinvention. Referring now to FIG. 5, a computer system 51 forimplementing the present invention can comprise, inter alia, a centralprocessing unit (CPU) 52, a memory 53 and an input/output (I/O)interface 54. The computer system 51 is generally coupled through theI/O interface 54 to a display 55 and various input devices 56 such as amouse and a keyboard. The support circuits can include circuits such ascache, power supplies, clock circuits, and a communication bus. Thememory 53 can include random access memory (RAM), read only memory(ROM), disk drive, tape drive, etc., or a combinations thereof. Thepresent invention can be implemented as a routine 57 that is stored inmemory 53 and executed by the CPU 52 to process the signal from thesignal source 58. As such, the computer system 51 is a general purposecomputer system that becomes a specific purpose computer system whenexecuting the routine 57 of the present invention.

The computer system 51 also includes an operating system and microinstruction code. The various processes and functions described hereincan either be part of the micro instruction code or part of theapplication program (or combination thereof) which is executed via theoperating system. In addition, various other peripheral devices can beconnected to the computer platform such as an additional data storagedevice and a printing device.

It is to be further understood that, because some of the constituentsystem components and method steps depicted in the accompanying figurescan be implemented in software, the actual connections between thesystems components (or the process steps) may differ depending upon themanner in which the present invention is programmed. Given the teachingsof the present invention provided herein, one of ordinary skill in therelated art will be able to contemplate these and similarimplementations or configurations of the present invention.

While the present invention has been described in detail with referenceto exemplary embodiments, those skilled in the art will appreciate thatvarious modifications and substitutions can be made thereto withoutdeparting from the spirit and scope of the invention as set forth in theappended claims.

APPENDIX Nomenclature

The main notation used throughout the present paper is summarized below.

Two-Letter Mnemonic of Energy Resources

BL Building.

BR Branch.

CD Condenser.

CG Chiller group.

CH Chiller.

CM Building group

DA Deaeretor.

EG Electric grid (excluding buildings)

GN Generator.

GS Gas

GT Gas turbine.

HR Heat recovery steam generator

IF Building-grid interface

LD Load

ND Node.

PM Pump.

ST Steam turbine.

TE Terminal.

TS Thermal energy storage (TES).

WR Water

Parameters

-   -   π_(igk) ^(CG) Chilled water supply temperature of combined        stream from all chillers (K)    -   π_(1ck) ^(CD) Condenser i pressure 1 in Steam Loop (kPa)    -   π_(2ck) ^(CD) Condenser i pressure 2 in Steam Loop (kPa)    -   π_(1ck) ^(CH) Power consumed by chiller c condenser pumps (kW)    -   π_(2ck) ^(CH) Power consumed by chiller c chilled water pumps        (kW)    -   π_(3ck) ^(CH) Power requirement by chiller c cooling tower fans        (kW)    -   π_(1mk) ^(CM) Cost of electricity purchase from grid by building        group m at time k ($/KWh)    -   π_(2mk) ^(CM) Electricity demand from building group m at time k        (MW)    -   π_(3mk) ^(CM) Cooling demand from building group m at time k        (kW)    -   π_(1dk) ^(DA) Deaerator d pressure in Steam Loop (kPa)    -   π_(1gk) ^(GS) Ideal gas specific heat capacity of steam        (KJ/kg-K)    -   π_(1gk) ^(CT) Cost of fuel purchase for GT i at time k ($/kg)    -   π_(1hk) ^(HR) Water stream pressure loss factor in HR h at time        k    -   π_(2hk) ^(HR) Overall effectiveness of HR h at time k    -   π_(1sk) ^(ST) Overall efficiency of ST s at time k    -   π_(1tk) ^(TS) Thermal transport delay time in discharging mode        of TS t at time k(s)    -   π_(2tk) ^(TS) TES thermal transport delay time in charging        mode(s)    -   π_(3tk) ^(TS) Time constant multiplier associated with TS upper        layer in charging mode.    -   π_(4tk) ^(TS) Time constant multiplier associated with the TS        upper layer in discharging mode.    -   π_(5tk) ^(TS) Inter-layer heat transfer coefficient of TS in        discharging mode (kW/m²-K)    -   π_(6tk) ^(TS) Time constant multiplier associated with the TES        bottom layer in discharging mode.    -   π_(7tk) ^(TS) Inter-layer heat transfer coefficient of TES in        charging mode (kW/m²-K)    -   π_(8tk) ^(TS) Time constant multiplier associated with the TES        bottom layer in charging mode.    -   π_(9tk) ^(TS) Area of TES tank (m²)    -   π_(1wk) ^(WR) Density of water (kg/s)    -   π_(2wk) ^(WR) Enthalpy of vaporization of water.    -   π_(2wk) ^(WR) Specific heat capacity of water (kJ/kg-K)        Decision Variables    -   x_(1gk) ^(CG) Temperature of return water to chiller group g at        time k (K)    -   x_(2gk) ^(CG) Net chilled water mass flow rate through chiller        group g at time k (kg/s)    -   x_(1ck) ^(CH) Power requirement by chiller compressor c at time        k (kW)    -   x_(2ck) ^(CH) Supply water temperature from chiller c at time k        (K)    -   x_(3ck) ^(CH) Cooling provided by chiller c at time k (kW)    -   x_(4ck) ^(CH) Chilled water mass flow rate in chiller c at time        k (kg/s)    -   x_(1mk) ^(CM) Chilled water flow rate through building group m        at time k (kg/s)    -   x_(2mk) ^(CM) Chilled water supply temperature to building group        m at time k (K)    -   x_(3mk) ^(CM) Temperature of return water from building group m        at time k (K)    -   x_(4ck) ^(CM) Mass flow rate of steam used for heating building        group m at time k (kg/s)    -   x_(5mk) ^(CM) Power purchase from electric grid by building        group m at time k (MW)    -   x_(1gk) ^(GT) Flow rate of exhaust gas of GT g at time k (kg/s)    -   x_(2gk) ^(GT) Turbine exit temperature from GT g at time k (K)    -   x_(3gk) ^(GT) Power produced by GT g at time k (MW)    -   x_(4gk) ^(GT) Fuel consumption rate by GT g at time k (kg/s)    -   x_(1hk) ^(HR) Inlet water temperature to HR h at time k (K)    -   x_(2hk) ^(HR) Steam temperature exiting HR h at time k (K)    -   x_(3hk) ^(HR) Saturation temperature in HR h at time k.    -   x_(4hk) ^(HR) Gas temperature exiting HR h at time k (K)    -   x_(5hk) ^(HR) Steam pressure exiting HR h at time k (kPa)    -   x_(6hk) ^(HR) Saturation pressure in HR h at time k.    -   x_(7hk) ^(HR) Steam pressure entering HR h at time k (kPa).    -   x_(8hk) ^(HR) Water mass flow rate through HR h at time k (kg/s)    -   x_(1pk) ^(PM) Power consumption by pump p in steam loop at time        k (kW)    -   x_(1sk) ^(ST) Power produced by ST s at time k (MW)    -   x_(2sk) ^(ST) Mass flow rate of steam used to drive ST s at time        k (kg/s)    -   x_(2sk) ^(ST) Fraction of cogenerated steam used to drive ST s        at time k    -   x_(1tk) ^(TS) Top layer temperature in 2-zone model of TS t at        time k (K)    -   x_(2tk) ^(TS) Bottom layer temperature in 2-zone model of TS t        at time k (K)    -   x_(3tk) ^(TS) Inlet stream temperature in charging mode of TS t        at time k (K)    -   x_(4tk) ^(TS) Outlet stream temperature in charging mode of TS t        at time k (K)    -   x_(5tk) ^(TS) Outlet stream temperature in discharging mode of        TS t at time k (K)    -   x_(6tk) ^(TS) Inlet stream temperature in discharging mode of TS        t at time k (K)    -   x_(7tk) ^(TS) Chilled water circulation through TS t at time k        (kg/s)    -   x_(8tk) ^(TS) Outlet stream temperature in charging mode of TS t        at time k, as predicted by detailed model (K)    -   x_(9tk) ^(TS) Outlet stream temperature in discharging mode of        TS t at time k, as predicted by detailed model (K)

What is claimed is:
 1. A computer-implemented method for optimizing acost of electric power generation in a smart site energy managementmodel, the method implemented by a computer comprising the steps of:providing, to a computer system, a cost function ζ that models a smartbuilding-grid energy system of a plurality of buildings on a siteinterconnected with electric power grid energy resources and constraintsdue to a building model, an electric grid model, and a building-gridinterface model, wherein decision variables for each of the buildingmodel, the electric grid model, and the building-grid interface modelare box-constrained; and minimizing, by the computer system, said costfunction subject to the building model constraints, the electric gridmodel constraints, and building-grid interface model constraints,wherein the building model includes a plurality of buildings, one ormore gas turbines and generators as on-site sources for electricitypower, electric chillers, pumps and thermal energy storage tanks,wherein waste heat in exhaust gas from the gas turbines is used togenerate steam in one or more heat recovery steam generator (HRSG)units, steam from the HRSG units drives one or more absorption chillersto generate cooling, one or more steam turbines to drive additionalelectricity generators, or provides heating needs for the plurality ofbuildings, the electric grid includes one or more loads and is modeledas an undirected graph G=(N, E), where N and E denote sets of electricnodes and branches, respectively, wherein each branch has an originterminal and a destination terminal, each terminal is connected to anelectric node, a voltage at node n₁ is a complex number defined by itsmagnitude and angle, and a power flow at terminal t₁ is a complex numberdefined by its real and imaginary parts, and the building-grid interfacemodel constrains a building cooling load to be satisfied by chilleroperations at all times, and an electricity load to be satisfied byon-site generation and the electric grid.
 2. The method of claim 1,wherein the cost function is${\varsigma = {\sum\limits_{k}{\sum\limits_{j}\left\lbrack {\rho_{jk}^{T_{j}}{u_{djk}^{T_{j}}\left( x_{jk}^{T_{j}} \right)}} \right\rbrack^{2}}}},$wherein ρ is a vector of parameters of components of the building modeland electric grid model, u_(d) is a basis vector of d-degreepolynomials, x is a vector of decision variables of the building modeland electric grid model, wherein a dimension of ρ and u_(d)(x) is givenby ${{u_{d}} = \begin{pmatrix}{{x} + d} \\d\end{pmatrix}},$ j is an index over the components of the building modeland electric grid model, T_(j) denotes one of the components of thebuilding model and electric grid model, and k is an index of timemeasurements of said parameters and decision variable values, and saidconstraints are P^(BL)(x)=0, P^(EG)(x)=0, P^(IF)(x)=0, and x≧0, whereinP^(BL)(x), P^(EG)(x), and P^(IF)(x) are sets of polynomials defining thebuilding (BL) model constraints, the electric grid (EG) modelconstraints, and the building-grid interface model (IF) constraints,respectively.
 3. The method of claim 2, wherein the components of thebuilding model include the electric and absorption chillers, a group ofchillers, a group of the plurality of buildings, the generators, the gasturbines, the HRSG units, the pumps, the steam turbines, the thermalenergy storages, and the components of the electric grid model includeloads, nodes, and terminals.
 4. The method of claim 2, wherein said costfunction is reformulated in a semi-definite programming format as${\min\limits_{X}{C \cdot X}},{{s.t.\mspace{14mu}{A^{BL}(X)}} = b^{BL}}$A^(EG)(X) = b^(EG) A^(IF)(X) = b^(IF) X ≽ 0, rank(X) = 1, byreformulating non-linear terms in the cost function and constraints asbilinear terms and using a matrix inner product operator on the bilinearterms, wherein C is a matrix formed of the parameters of the componentsof the building model and electric grid model, X is a matrix formed ofthe decision variables of the building model and electric grid model,A^(BL)(x), A^(EG)(X), and A^(IF)(X) are reformulations of the building(BL) model constraints, the electric grid (EG) constraints, and thebuilding-grid interface model (IF) constraints, respectively, in asemi-definite programming format, b^(BL), b^(EG), and b^(IF) areconstants derived from the constraints of the building model andelectric grid model, and $X = {\begin{bmatrix}1 & x^{T} \\x & {xx}^{T}\end{bmatrix}.}$
 5. The method of claim 4, further comprising relaxingthe condition rank(X)=1, and solving a following primal semi-definiteprogram${{\min\limits_{X}{C \cdot X}} - {\mu\;{\ln\left( {\det(X)} \right)}}},{{wherein}\mspace{14mu}\mu\mspace{14mu}{is}\mspace{14mu} a\mspace{14mu}{Lagrange}\mspace{14mu}{multiplier}},{{{subject}\mspace{14mu}{to}\mspace{14mu}{A^{BL}(X)}} = b^{BL}},{{A^{EG}(X)} = b^{EG}},{{A^{IF}(X)} = b^{IF}},{X \succ 0},$for positive real values of μ as μ approaches zero.
 6. The method ofclaim 5, wherein solving the primal semi-definite program comprises:imposing optimality conditions on the primal semi-definite program toderiveC−Z−

^(BL)(y ^(BL))−

^(EG)(y ^(EG))−

^(IF)(y ^(IF))=0,

^(BL)(X)−b ^(BL)=0,

^(EG)(X)−b ^(EG)=0,

^(IF)(X)−b ^(IF)=0,XZ−μI=0, wherein Z=μX⁻¹, X

0, and Ã^(BL)(y^(BL)), Ã^(EG)(y^(EG)) and Ã^(IF)(y^(IF)) are transposesof A^(BL)(y^(BL)), A^(EG)(y^(EG)) and A^(IF)(y^(IF)) where y^(BL),y^(EQ) and y^(IF) are variables in a dual program of the primalsemi-definite program; initializing a solution (X, y, Z) to an initialpoint (X⁰, y⁰, Z⁰) wherein X⁰

0, Z⁰

0, and y is a vector that aggregates y^(BL), y^(EQ), and y^(IF);choosing a search direction (ΔX, Δy, ΔZ), wherein ΔX, Δy, ΔZ representschanges in X, y, Z; choosing a primal step length α_(p) and a dual steplength α_(d) that satisfy X+α_(p)ΔX

0 and Z+α_(d)ΔZ

0; and updating X←X+α_(p)ΔX and (y, Z)←α_(d)(Δy, ΔZ).
 7. The method ofclaim 6, further comprising repeating the steps of choosing a searchdirection, choosing a primal step length and a dual step length, andupdating a current iterate (X, y, Z) until the current iterate satisfiesa stopping condition.
 8. The method of claim 6, wherein choosing asearch direction comprises solvingΔZ+

^(BL)(Δy ^(BL))+

^(EG)(Δy ^(EG))+

^(IF)(Δy ^(IF))=−r ^(X),

^(BL)(ΔX)=−r ^(BL),

^(EG)(ΔX)=−r ^(EG),

^(IF)(ΔX)=−r ^(IF),

_(P)((ΔX)Z+XΔZ)=−r ^(XZ),whereinr ^(X) =C−ΔZ−

^(BL)(Δy ^(BL))−

^(EG)(Δy ^(EG))−

^(IF)(Δy ^(IF)),r ^(BL)=

^(BL)(X)−b ^(BL),r ^(EG)=

^(EG)(X)−b ^(EG),r ^(IF)=

^(IF)(X)−b ^(IF),r ^(XZ)=

_(P)(XZ)−σμI, wherein I is an identity matrix,σε[0,1],H _(P)(XZ)=½(PXZP ⁻¹+(PXZP ⁻¹)^(T)), andH _(P)(ΔXZ+XΔZ)=½(P(ΔXZ+XΔZ)P ⁻¹+(P(ΔXZ+XΔZ)P ⁻¹)^(T)) for anon-singular matrix P.
 9. A computer-implemented method for optimizing acost of electric power generation in a smart site energy managementmodel, the method implemented by a computer system comprising the stepsof: providing, to the computer system, a semidefinite program${\min\limits_{X}{C \cdot X}} - {\mu\;{\ln\left( {\det(X)} \right)}}$subject  to  A^(BL)(X) = b^(BL), A^(EG)(X) = b^(EG), A^(IF)(X) = b^(IF), X ≻ 0,wherein C·X is a cost function that models a smart building-grid energysystem of a plurality of buildings on a site interconnected withelectric power grid energy resources wherein C is a matrix formed ofparameters of components of a building model and an electric grid model,X is a matrix formed of decision variables of the building model andelectric grid model, μ is a Lagrange multiplier, A^(BL)(X), A^(EG)(X),and A^(IF)(X) are constraints due to a building (BL) model, an electricgrid (EG) model, and a building-grid interface model (IF), respectively,in a semi-definite programming format, b^(BL), b^(EG), and b^(IF) areconstants derived from the constraints of the building model andelectric grid model, and ${X = \begin{bmatrix}1 & x^{T} \\x & {xx}^{T}\end{bmatrix}},$ wherein a rank condition on X has been relaxed; andsolving, by the computer system,${{\min\limits_{X}{C \cdot X}} - {\mu\;{\ln\left( {\det(X)} \right)}}},$subject to the constraints, for positive real values of μ as μapproaches zero.
 10. The method of claim 9, wherein C·X is derived byreformulating non-linear terms in a function ζ and constraintsP^(BL)(x)=0, P^(EG)(x)=0, P^(IF)(x)=0, and x≧0, as bilinear terms, usinga matrix inner product operator on the bilinear terms to reformulate thebilinear terms in a semi-definite programming format, wherein${\varsigma = {\sum\limits_{k}{\sum\limits_{j}\left\lbrack {\rho_{jk}^{T_{j}}{u_{djk}^{T_{j}}\left( x_{jk}^{T_{j}} \right)}} \right\rbrack^{2}}}},$wherein ρ is a vector of parameters of components of the building modeland electric grid model, u_(d) is a basis vector of d-degreepolynomials, x is a vector of decision variables of the building modeland electric grid model, wherein a dimension of ρ and u_(d)(x) is givenby ${{u_{d}} = \begin{pmatrix}{{x} + d} \\d\end{pmatrix}},$ j is an index over the components of the building modeland electric grid model, T_(j) denotes one of the components of thebuilding model and electric grid model, and k is an index of timemeasurements of said parameters and decision variable values, andP^(BL)(x), P^(EG)(x), and P^(IF)(x) are sets of polynomials defining thebuilding (BL) model constraints, the electric grid (EG) modelconstraints, and the building-grid interface model (IF) constraints,respectively.
 11. The method of claim 9, wherein decision variables foreach of the building model, the electric grid model, and thebuilding-grid interface model are box-constrained.
 12. The method ofclaim 9, wherein the building model includes a plurality of buildings,one or more gas turbines and generators as on-site sources forelectricity power, electric chillers, pumps and thermal energy storagetanks, wherein waste heat in exhaust gas from the gas turbines is usedto generate steam in one or more heat recovery steam generator (HRSG)units, steam from the HRSG units drives one or more absorption chillersto generate cooling, steam turbines to drive additional electricitygenerators, or provides heating needs for the plurality of buildings,the electric grid includes one or more loads and is modeled as anundirected graph G=(N, E), where N and E denote sets of electric nodesand branches, respectively, wherein each branch has an origin terminaland a destination terminal, each terminal is connected to an electricnode, a voltage at node n₁ is a complex number defined by its magnitudeand angle, and a power flow at terminal t₁ is a complex number definedby its real and imaginary parts, and the building-grid interface modelconstrains a building cooling load to be satisfied by chiller operationsat all times, and an electricity load to be satisfied by on-sitegeneration and the electric grid.
 13. The method of claim 9, whereinsolving the primal semi-definite program comprises: imposing optimalityconditions on the primal semi-definite program to deriveC−Z−

^(BL)(y ^(BL))−

^(EG)(y ^(EG))−

^(IF)(y ^(IF))=0,

^(BL)(X)−b ^(BL)=0,

^(EG)(X)−b ^(EG)=0,

^(IF)(X)−b ^(IF)=0,XZ−μI=0, wherein Z=μX⁻¹, X

0, and Ã^(BL)(y^(BL)), Ã^(EG)(y^(EG)) and Ã^(IF)(y^(IF)) are transposesof A^(BL)(y^(BL)), A^(EG)(y^(EG)) and A^(IF)(y^(IF)) where y^(BL),y^(EQ) and y^(IF) are variables in a dual program of the primalsemi-definite program; initializing a solution (X, y, Z) to an initialpoint (X⁰, y⁰, Z⁰) wherein X⁰

0, Z⁰

0, and y is a vector that aggregates y^(BL), y^(EQ), and y^(IF);choosing a search direction (ΔX, Δy, ΔZ), wherein ΔX, Δy, ΔZ representschanges in X, y, Z; choosing a primal step length α_(p) and a dual steplength α_(d) that satisfy X+α_(p)ΔX

0 and Z+α_(d)ΔZ

0; and updating X←X+α_(p)ΔX and (y, Z)←α_(d)(Δy, ΔZ); and repeating thesteps of choosing a search direction, choosing a primal step length anda dual step length, and updating a current iterate (X, y, Z) until thecurrent iterate satisfies a stopping condition.
 14. The method of claim13, wherein choosing a search direction comprises solvingΔZ+

^(BL)(Δy ^(BL))+

^(EG)(Δy ^(EG))+

^(IF)(Δy ^(IF))=−r ^(X),

^(BL)(ΔX)=−r ^(BL),

^(EG)(ΔX)=−r ^(EG),

^(IF)(ΔX)=−r ^(IF),

_(P)((ΔX)Z+XΔZ)=−r ^(XZ),whereinr ^(X) =C−ΔZ−

^(BL)(Δy ^(BL))−

^(EG)(Δy ^(EG))−

^(IF)(Δy ^(IF)),r ^(BL)=

^(BL)(X)−b ^(BL),r ^(EG)=

^(EG)(X)−b ^(EG),r ^(IF)=

^(IF)(X)−b ^(IF),r ^(XZ)=

_(P)(XZ)−σμI, wherein I is an identity matrix,σε[0,1],H _(P)(XZ)=½(PXZP ⁻¹+(PXZP ⁻¹)^(T)), andH _(P)(ΔXZ+XΔZ)=½(P(ΔXZ+XΔZ)P ⁻¹+(P(ΔXZ+XΔZ)P ⁻¹)^(T)) for anon-singular matrix P.
 15. A non-transitory program storage devicereadable by a computer, tangibly embodying a program of instructionsexecuted by the computer to perform the method steps for optimizing acost of electric power generation in a smart site energy managementmodel, comprising the steps of: providing a cost function ζ that modelsa smart building-grid energy system of a plurality of buildings on asite interconnected with electric power grid energy resources andconstraints due to a building model, an electric grid model, and abuilding-grid interface model, wherein decision variables for each ofthe building model, the electric grid model, and the building-gridinterface model are box-constrained; and minimizing said cost functionsubject to the building model constraints, the electric grid modelconstraints, and building-grid interface model constraints, wherein thebuilding model includes a plurality of buildings, one or more gasturbines and generators as on-site sources for electricity power,electric chillers, pumps and thermal energy storage tanks, wherein wasteheat in exhaust gas from the gas turbines is used to generate steam inone or more heat recovery steam generator (HRSG) units, steam from theHRSG units drives one or more absorption chillers to generate cooling,one or more steam turbines to drive additional electricity generators,or provides heating needs for the plurality of buildings, the electricgrid includes one or more loads and is modeled as an undirected graphG=(N, E), where N and E denote sets of electric nodes and branches,respectively, wherein each branch has an origin terminal and adestination terminal, each terminal is connected to an electric node, avoltage at node n₁ is a complex number defined by its magnitude andangle, and a power flow at terminal t₁ is a complex number defined byits real and imaginary parts, and the building-grid interface modelconstrains a building cooling load to be satisfied by chiller operationsat all times, and an electricity load to be satisfied by on-sitegeneration and the electric grid.
 16. The computer readable programstorage device of claim 15, wherein the cost function is${\varsigma = {\sum\limits_{k}{\sum\limits_{j}\left\lbrack {\rho_{jk}^{T_{j}}{u_{djk}^{T_{j}}\left( x_{jk}^{T_{j}} \right)}} \right\rbrack^{2}}}},$wherein ρ is a vector of parameters of components of the building modeland electric grid model, u_(d) is a basis vector of d-degreepolynomials, x is a vector of decision variables of the building modeland electric grid model, wherein a dimension of ρ and u_(d)(x) is givenby ${{u_{d}} = \begin{pmatrix}{{x} + d} \\d\end{pmatrix}},$ j is an index over the components of the building modeland electric grid model, T_(j) denotes one of the components of thebuilding model and electric grid model, and k is an index of timemeasurements of said parameters and decision variable values, and saidconstraints are P^(BL)(x)=0, P^(EG)(x)=0, P^(IF)(x)=0, and x≧0, whereinP^(BL)(x), P^(EG)(x), and P^(IF)(x) are sets of polynomials defining thebuilding (BL) model constraints, the electric grid (EG) modelconstraints, and the building-grid interface model (IF) constraints,respectively.
 17. The computer readable program storage device of claim16, wherein the components of the building model include the electricand absorption chillers, a group of chillers, a group of the pluralityof buildings, the generators, the gas turbines, the HRSG units, thepumps, the steam turbines, the thermal energy storages, and thecomponents of the electric grid model include loads, nodes, andterminals.
 18. The computer readable program storage device of claim 16,wherein said cost function is reformulated in a semi-definiteprogramming format as${\min\limits_{X}{C \cdot X}},{{s.t.\mspace{14mu}{A^{BL}(X)}} = b^{BL}}$A^(EG)(X) = b^(EG) A^(IF)(X) = b^(IF) X ≽ 0, rank(X) = 1, byreformulating non-linear terms in the cost function and constraints asbilinear terms and using a matrix inner product operator on the bilinearterms, wherein C is a matrix formed of the parameters of the componentsof the building model and electric grid model, X is a matrix formed ofthe decision variables of the building model and electric grid model,A^(BL)(X), A^(EG)(X), and A^(IF)(X) are reformulations of the building(BL) model constraints, the electric grid (EG) constraints, and thebuilding-grid interface model (IF) constraints, respectively, in asemi-definite programming format, b^(BL), b^(EG), and b^(IF) areconstants derived from the constraints of the building model andelectric grid model, and $X = {\begin{bmatrix}1 & x^{T} \\x & {xx}^{T}\end{bmatrix}.}$
 19. The computer readable program storage device ofclaim 18, the method further comprising relaxing the conditionrank(X)=1, and solving a following primal semi-definite program${{\min\limits_{X}{C \cdot X}} - {\mu\;{\ln\left( {\det(X)} \right)}}},\;{{wherein}\mspace{14mu}\mu\mspace{14mu}{is}\mspace{14mu} a\mspace{14mu}{Lagrange}\mspace{14mu}{multiplier}},{{{subject}\mspace{14mu}{to}\mspace{14mu}{A^{BL}(X)}} = b^{BL}},{{A^{EG}(X)} = b^{EG}},{{A^{IF}(X)} = b^{IF}},{X \succ 0},$for positive real values of μ as μ approaches zero.
 20. The computerreadable program storage device of claim 19, wherein solving the primalsemi-definite program comprises: imposing optimality conditions on theprimal semi-definite program to deriveC−Z−

^(BL)(y ^(BL))−

^(EG)(y ^(EG))−

^(IF)(y ^(IF))=0,

^(BL)(X)−b ^(BL)=0,

^(EG)(X)−b ^(EG)=0,

^(IF)(X)−b ^(IF)=0,XZ−μI=0, wherein Z=μX⁻¹, X

0, and Ã^(BL)(y^(BL)), Ã^(EG)(y^(EG)) and Ã^(IF)(y^(IF)) are transposesof A^(BL)(y^(BL)), A^(EG)(y^(EG)) and A^(IF)(y^(IF)) where y^(BL),y^(EQ) and y^(IF) are variables in a dual program of the primalsemi-definite program; initializing a solution (X, y, Z) to an initialpoint (X⁰, y⁰, Z⁰) wherein X⁰

0, Z⁰

0, and y is a vector that aggregates y^(BL), y^(EQ), and y^(IF);choosing a search direction (ΔX, Δy, ΔZ), wherein ΔX, Δy, ΔZ representschanges in X, y, Z; choosing a primal step length α_(p) and a dual steplength α_(d) that satisfy X+α_(p)ΔX

0 and Z+α_(d)ΔZ

0; and updating X←X+α_(p)ΔX and (y, Z)←α_(d)(Δy, ΔZ).
 21. The computerreadable program storage device of claim 20, the method furthercomprising repeating the steps of choosing a search direction, choosinga primal step length and a dual step length, and updating a currentiterate (X, y, Z) until the current iterate satisfies a stoppingcondition.
 22. The computer readable program storage device of claim 20,wherein choosing a search direction comprises solvingΔZ+

^(BL)(Δy ^(BL))+

^(EG)(Δy ^(EG))+

^(IF)(Δy ^(IF))=−r ^(X),

^(BL)(ΔX)=−r ^(BL),

^(EG)(ΔX)=−r ^(EG),

^(IF)(ΔX)=−r ^(IF),

_(P)((ΔX)Z+XΔZ)=−r ^(XZ),whereinr ^(X) =C−ΔZ−

^(BL)(Δy ^(BL))−

^(EG)(Δy ^(EG))−

^(IF)(Δy ^(IF)),r ^(BL)=

^(BL)(X)−b ^(BL),r ^(EG)=

^(EG)(X)−b ^(EG),r ^(IF)=

^(IF)(X)−b ^(IF),r ^(XZ)=

_(P)(XZ)−σμI, wherein I is an identity matrix,σε[0,1],H _(P)(XZ)=½(PXZP ⁻¹+(PXZP ⁻¹)^(T)), andH _(P)(ΔXZ+XΔZ)=½(P(ΔXZ+XΔZ)P ⁻¹+(P(ΔXZ+XΔZ)P ⁻¹)^(T)) for anon-singular matrix P.